library(parallel)
library(doMC)
library(DHARMa)
library(broom.mixed)
library(gratia)
library(MuMIn)
library(mgcv)
library(bbmle)
library(lme4)
library(merTools)
library(magrittr)
library(gridExtra)
library(ggplot2); theme_set(theme_classic())
library(tidymv)
library(purrr)
library(plyr)
library(dplyr)

1 - Variables

table 1: Variáveis utilizadas na descrição estatística

Variables
code description class range
n_nRef number of SADs not refuted by the KS bootstrap test with critical p-value 0.05 continuous (0 ; 100)
diffS_mean diffS = (S_MN - S_obs) / S_obs, diffS_mean = mean(diffS); S_MN = species richness estimated by MN continuous (-0.977 ; 4.78)
U_med average of 10 replicates of the speciation rate estimated by MNEE continuous (8.860e-5 ; 4.221e-2)
p proportion of tree cover continuous (0.013 ; 0.961)
MN Neutral Model (EE = spatial explicit; EI = spatial implict) category 2 levels
d mean dispersal distance (meters) continuous (0.456 ; 19.183)
k proportion of propagules until the nearest neighborhood seq(0.99 : 0.05) category 20 levels
d_Lplot mean dispersal distance / side of the sample area continuous (0.02 ; 0.192)
S_obs observed species richness integer (26 ; 458)
Ntotal number of individuals in the sample area integer (649 ; 12105)
SiteCode sample area code category 103 levels

Figure 1 Possible Predictor Variables

2 - Number of unrefuted SADs

Figura 2.1 número de SADs não refutadas ~ p * k * MN. A linha azul é uma estimativa baseada em ‘loess’.

Figura 2.2 número de SADs não refutadas ~ (d * MN|Site ~ p). Site está ordenado pelo valor de p. 

2.1 GLMM binomial

Modelo cheio 1:

linear term

  1. ~ (p + p^2) * (d + d^2) * MN; if : var_dispersao = contiguous
  2. ~ ( p + p^2) * k * MN; if : var_dispersao = category

random term

  1. 1|SiteCode
  2. MN|SiteCode
  3. (var_dispersao + var_dispersao^2)*MN|SiteCode, if : var = dispersao : contiguous
# dados
df_md <- df_resultados %>% 
  select(n_nRef,p,SiteCode,MN,p.z,k,d.z,d_Lplot.z)
df_md[,6:8] <- apply(df_md[,6:8],2,as.character)
df_md %<>% 
  gather(key = var_dispersao, value = value_dispersao, k:d_Lplot.z) 
# funcao ajuste dos modelos
f_md <- function(dados){
  var_dispersao <- unique(dados$var_dispersao)
  if(var_dispersao == "k"){
    l_md <- vector("list",2)
    names(l_md) <- c(paste0(var_dispersao," 1|Site"),
                     paste0(var_dispersao," MN|Site"))
    dados$value_dispersao <- factor(dados$value_dispersao,levels = unique(dados$value_dispersao)[20:1])
    l_md[[1]] <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                         (p.z + I(p.z^2)) * value_dispersao * MN +
                         (1|SiteCode),
                       family = "binomial",data=dados,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
    l_md[[2]] <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                         (p.z + I(p.z^2)) * value_dispersao * MN +
                         (MN|SiteCode),
                       family = "binomial",data=dados,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
  }else{
    dados$value_dispersao <- as.numeric(dados$value_dispersao)
    l_md <- vector("list",3)
    names(l_md) <- c(paste0(var_dispersao," 1|Site"),
                     paste0(var_dispersao," MN|Site"),
                     paste0(var_dispersao," (",var_dispersao,"+",var_dispersao,"^2)*MN|Site"))
    l_md[[1]] <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                         (p.z + I(p.z^2)) * (value_dispersao + I(value_dispersao^2)) * MN + 
                         (1|SiteCode),
                       family = "binomial",data=dados,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
    l_md[[2]] <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                         (p.z + I(p.z^2)) * (value_dispersao + I(value_dispersao^2)) * MN + 
                         (MN|SiteCode),
                       family = "binomial",data=dados,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
    l_md[[3]] <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                         (p.z + I(p.z^2)) * (value_dispersao + I(value_dispersao^2)) * MN + 
                         ( (value_dispersao + I(value_dispersao^2) ) * MN|SiteCode),
                       family = "binomial",data=dados,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
    
  }
  return(l_md)
}
registerDoMC(3)
l_nRef__modeloCheio <- dlply(df_md,.(var_dispersao),f_md,.parallel = TRUE)
names(l_nRef__modeloCheio) <- NULL
l_nRef__modeloCheio <- do.call(c,l_nRef__modeloCheio)
#
f_warning <- function(glmm_object){
  # glmm_object <- l_nRef__modeloCheio[[4]]
  v_message <- glmm_object@optinfo$conv$lme4$messages %>% 
    as.character()
  if(length(v_message)==0){v_message <- "OK"}
  if(length(v_message)>1){v_message <- v_message[1] }
  return(v_message)
}
df_auditoria.md <- ldply(l_nRef__modeloCheio,f_warning,.id="glmer") %>% 
  rename(warning_message = V1)

Table 2.1 Modelos Cheios estimados e avisos de convergência

glmer warning_message
d_Lplot.z 1|Site Model failed to converge with max|grad| = 0.00343548 (tol = 0.002, component 1)
d_Lplot.z MN|Site Model failed to converge with max|grad| = 0.00836802 (tol = 0.002, component 1)
d_Lplot.z (d_Lplot.z+d_Lplot.z^2)*MN|Site OK
d.z 1|Site Model failed to converge with max|grad| = 40.1827 (tol = 0.002, component 1)
d.z MN|Site Model failed to converge with max|grad| = 0.00746139 (tol = 0.002, component 1)
d.z (d.z+d.z^2)*MN|Site OK
k 1|Site OK
k MN|Site OK

allFit

i <- 1
names(l_nRef__modeloCheio[v_glmerUpdate])[i]
#
## update1:
md_allFit <- allFit(l_nRef__modeloCheio[v_glmerUpdate][[i]],maxfun = 1e5, parallel = 'multicore', ncpus = 3)

1) d_Lplot 1|Site

7 optimizer(s) failed

2) d_Lplot MN|Site

7 optimizer(s) failed

3) d 1|Site

7 optimizer(s) failed

4) d MN|Site

7 optimizer(s) failed

Modelo cheio 2:

# dados
df_md <- df_resultados %>% 
  select(n_nRef,p,SiteCode,MN,p.z,k,d.z,d_Lplot.z)
df_md[,6:8] <- apply(df_md[,6:8],2,as.character)
df_md %<>% 
  gather(key = var_dispersao, value = value_dispersao, k:d_Lplot.z) 
# funcao ajuste dos modelos
f_md <- function(dados){
  var_dispersao <- unique(dados$var_dispersao)
  if(var_dispersao == "k"){
    l_md <- vector("list",2)
    names(l_md) <- c(paste0(var_dispersao," 1|Site"),
                     paste0(var_dispersao," MN|Site"))
    dados$value_dispersao <- factor(dados$value_dispersao,levels = unique(dados$value_dispersao)[20:1])
    l_md[[1]] <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                         (p.z + I(p.z^2)) * value_dispersao * MN +
                         (1|SiteCode),
                       family = "binomial",data=dados,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
    l_md[[2]] <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                         (p.z + I(p.z^2)) * value_dispersao * MN +
                         (MN|SiteCode),
                       family = "binomial",data=dados,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
  }else{
    dados$value_dispersao <- as.numeric(dados$value_dispersao)
    l_md <- vector("list",4)
    names(l_md) <- c(paste0(var_dispersao," 1|Site"),
                     paste0(var_dispersao," MN|Site"),
                     paste0(var_dispersao," ",var_dispersao,"*MN|Site"),
                     paste0(var_dispersao,"^2 (",var_dispersao,"+",var_dispersao,"^2)*MN|Site"))
    l_md[[1]] <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                         (p.z + I(p.z^2)) * value_dispersao * MN + 
                         (1|SiteCode),
                       family = "binomial",data=dados,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
    l_md[[2]] <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                         (p.z + I(p.z^2)) * value_dispersao * MN + 
                         (MN|SiteCode),
                       family = "binomial",data=dados,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
    l_md[[3]] <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                         (p.z + I(p.z^2)) * value_dispersao * MN + 
                         (value_dispersao * MN|SiteCode),
                       family = "binomial",data=dados,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
    l_md[[4]] <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                         (p.z + I(p.z^2)) * (value_dispersao + I(value_dispersao^2)) * MN + 
                         ( (value_dispersao + I(value_dispersao^2) ) * MN|SiteCode),
                       family = "binomial",data=dados,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
    
  }
  return(l_md)
}
registerDoMC(3)
l_nRef__modeloCheio2 <- dlply(df_md,.(var_dispersao),f_md,.parallel = TRUE)
names(l_nRef__modeloCheio2) <- NULL
l_nRef__modeloCheio2 <- do.call(c,l_nRef__modeloCheio2)
#
f_warning <- function(glmm_object){
  v_message <- glmm_object@optinfo$conv$lme4$messages %>% 
    as.character()
  if(length(v_message)==0){v_message <- "OK"}
  if(length(v_message)>1){v_message <- v_message[1] }
  return(v_message)
}
df_auditoria.md <- ldply(l_nRef__modeloCheio2,f_warning,.id="glmer") %>% 
  rename(warning_message = V1)

Table 2.2 Modelos Cheios estimados e avisos de convergência

glmer warning_message
d_Lplot.z 1|Site Model failed to converge with max|grad| = 0.00284637 (tol = 0.002, component 1)
d_Lplot.z MN|Site OK
d_Lplot.z d_Lplot.z*MN|Site OK
d_Lplot.z^2 (d_Lplot.z+d_Lplot.z^2)*MN|Site OK
d.z 1|Site Model failed to converge with max|grad| = 0.00369235 (tol = 0.002, component 1)
d.z MN|Site OK
d.z d.z*MN|Site OK
d.z^2 (d.z+d.z^2)*MN|Site OK
k 1|Site OK
k MN|Site OK

allFit

d_Lplot.z 1|Site
7 optimizer(s) failed

d.z 1|Site
7 optimizer(s) failed

Comparação de Modelos Cheios:

Tabela 2.3 Comparação baseada em AICc dos modelos cheios

GLMM dAICc df weight
d^2 (d+d^2)*MN|Site 0.0 39 1
d_Lplot^2 (d_Lplot+d_Lplot^2)*MN|Site 274.6 39 <0.001
d d*MN|Site 38427.3 22 <0.001
d_Lplot d_Lplot*MN|Site 38586.5 22 <0.001
k MN|Site 74141.2 123 <0.001
d MN|Site 95710.0 15 <0.001
d_Lplot MN|Site 96279.6 15 <0.001
k 1|Site 106158.0 121 <0.001

Tabela 2.4 Coeficiente de Determinação Condicional e Marginal

##                   R2m       R2c
## theoretical 0.3232656 0.9128366
## delta       0.3111874 0.8787302

Diagnostico do modelo cheio mais plausível

Figura 2.3 Resíduos Quantílicos do modelo cheio plausível

Figura 2.4 Quantile-quantile plot random effects.

Figura 2.5 Predito e observado pelo modelo. Em vermelho linha y=x; em azul um modelo linear por sítio de amostragem.

Figura 2.6 Predito e observado para modelo cheio

Subset=MNEE

Estudo GAMM: smoother type

  1. Antes de comparar os modelos cheios preciso avaliar smoother type mais para utilizar na GAMM;
  2. Pelo observado nos dados (figuras 2.1 e 2.2) vou comparar dois smoother type: “tp” e “cr”
  3. Para isso vou ajuster um smoother para d por Sítio de amostragem com um mesmo parâmetro de penalização, dessa maneira eu considero a variação por sítio de amostragem na variável d sem um elevado custo computacional.
# dados
df_md <- df_resultados %>% filter(MN=="EE") %>% distinct()
df_md <- rbind(mutate(df_md,model_class = "cr"),
               mutate(df_md,model_class = "tp"))
df_md$model_class <- factor(df_md$model_class)
# função
f_md <- function(dados){
  model_class <- unique(dados$model_class)
  if(model_class == "cr"){
    md_ <- gam(cbind(n_nRef,100-n_nRef) ~ 
                 s(p.z,bs="cr") + s(d.z,bs="cr") + ti(p.z,d.z,bs=c("cr","cr")) +
                 s(d.z,SiteCode,bs="fs",xt=list(bs="cr"), m=1),
                data = dados, family = "binomial")
    
  }else{
    md_ <- gam(cbind(n_nRef,100-n_nRef) ~ 
                 s(p.z,bs="tp") + s(d.z,bs="tp") + ti(p.z,d.z,bs=c("tp","tp")) +
                 s(d.z,SiteCode,bs="fs",xt=list(bs="tp"), m=1),
                data = dados, family = "binomial")
    
  }
  return(md_)
}
registerDoMC(2)
l_nRefEE_estudoGAMM <- dlply(df_md,"model_class",f_md,.parallel = TRUE)

Comparação Modelos Cheios

# dados
df_md <- df_resultados %>% filter(MN=="EE") %>% 
  distinct() %>% 
  mutate(rep = 1)
df_md <- left_join(x=data.frame(rep=1,
                                md_class=c("glmm d+d^2|Site",
                                           "glmm d|Site",
                                           "glmm 1|Site",
                                           "gamm d|Site for each",
                                           "gamm d|Site common",
                                           "gamm 1|Site")),
                   y=df_md,
                   by="rep")
# função de ajuste
f_md <- function(dados){
  md_class <- unique(dados$md_class)
  if(md_class == "glmm d+d^2|Site"){
    md_ <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                      (p.z + I(p.z^2)) * (d.z + I(d.z^2)) + 
                      (d.z + I(d.z^2)|SiteCode),
                 family = "binomial",data=dados,
                 control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e6)),
                 na.action = "na.fail")
  }else if(md_class == "glmm d|Site"){
    md_ <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                      (p.z + I(p.z^2)) * (d.z + I(d.z^2)) + 
                      (d.z|SiteCode),
                 family = "binomial",data=dados,
                 control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e6)),
                 na.action = "na.fail")
  }else if(md_class == "glmm 1|Site"){
    md_ <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                      (p.z + I(p.z^2)) * (d.z + I(d.z^2)) + 
                      (1|SiteCode),
                 family = "binomial",data=dados,
                 control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e6)),
                 na.action = "na.fail")
  }else if(md_class == "gamm d|Site for each"){
    md_ <- gam(cbind(n_nRef,100-n_nRef) ~ 
                 s(p.z,bs="cr") + s(d.z,bs="cr") + ti(p.z,d.z,bs=c("cr","cr")) +
                 s(SiteCode,bs="re") + s(d.z,by=SiteCode,bs="cr",m=1),
                data = dados, family = "binomial")
  }else if(md_class == "gamm d|Site common"){
    md_ <- gam(cbind(n_nRef,100-n_nRef) ~ 
                 s(p.z,bs="cr") + s(d.z,bs="cr") + ti(p.z,d.z,bs=c("cr","cr")) +
                 s(d.z,SiteCode,bs="fs",xt=list(bs="cr"), m=1),
                data = dados, family = "binomial")
  }else{
    md_ <- gam(cbind(n_nRef,100-n_nRef) ~ 
                 s(p.z,bs="cr") + s(d.z,bs="cr") + ti(p.z,d.z,bs=c("cr","cr")) +
                 s(SiteCode,bs="re"),
                data = dados, family = "binomial")
  }
  return(md_)
}
registerDoMC(3)
l_nRefEE_mdCheio <- dlply(df_md,"md_class",f_md,.parallel = TRUE)
##                      dAICc   df               weight
## gamm d|Site for each     0.0 668.986466713628 1     
## gamm d|Site common     816.9 851.710244484507 <0.001
## glmm d+d^2|Site       5223.6 15               <0.001
## glmm d|Site           8188.0 12               <0.001
## gamm 1|Site          25233.8 126.840686107377 <0.001
## glmm 1|Site          27131.4 10               <0.001

Diagnostico modelo cheio

Figura 2.2.1 gam.check(md_nRefEE1)

## 
## Method: UBRE   Optimizer: outer newton
## full convergence after 33 iterations.
## Gradient range [-5.585756e-06,5.869852e-06]
## (score 0.6972787 & scale 1).
## eigenvalue range [-1.541502e-11,0.0008330883].
## Model rank =  1063 / 1065 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                               k'      edf k-index p-value
## s(p.z)                  9.00e+00 1.00e+00    1.34    1.00
## s(d.z)                  9.00e+00 8.92e+00    1.04    0.98
## ti(p.z,d.z)             1.60e+01 2.22e+00    0.99    0.56
## s(SiteCode)             1.03e+02 1.01e+02      NA      NA
## s(d.z):SiteCodeBAjuss   9.00e+00 3.52e+00    1.04    0.99
## s(d.z):SiteCodeBAlenc4  9.00e+00 6.21e+00    1.04    0.97
## s(d.z):SiteCodeBAuruc   9.00e+00 3.92e+00    1.04    0.98
## s(d.z):SiteCodeESsoor   9.00e+00 7.09e+00    1.04    0.97
## s(d.z):SiteCodeMGipia1  9.00e+00 5.60e+00    1.04    0.98
## s(d.z):SiteCodeMGlavr2  9.00e+00 7.39e+00    1.04    0.96
## s(d.z):SiteCodeMGlavr3  9.00e+00 1.00e+00    1.04    0.98
## s(d.z):SiteCodeMGlavr6  9.00e+00 1.77e+00    1.04    0.98
## s(d.z):SiteCodeMGubera  9.00e+00 8.90e+00    1.04    0.97
## s(d.z):SiteCodeMGuberl1 9.00e+00 3.39e+00    1.04    0.98
## s(d.z):SiteCodeMGuberl3 9.00e+00 2.48e+00    1.04    0.97
## s(d.z):SiteCodeMGuberl5 9.00e+00 8.91e+00    1.04    0.96
## s(d.z):SiteCodeMGuberl6 9.00e+00 5.25e+00    1.04    0.98
## s(d.z):SiteCodeMGuberl7 9.00e+00 4.38e+00    1.04    0.98
## s(d.z):SiteCodeMGvico1  9.00e+00 8.53e+00    1.04    0.97
## s(d.z):SiteCodeMGvico16 9.00e+00 8.15e+00    1.04    0.98
## s(d.z):SiteCodeMGvico3  9.00e+00 7.35e+00    1.04    0.98
## s(d.z):SiteCodeMGvico4  9.00e+00 8.13e+00    1.04    0.96
## s(d.z):SiteCodeMSdour   9.00e+00 1.63e-06    1.04    0.97
## s(d.z):SiteCodePEalia   9.00e+00 6.03e+00    1.04    0.98
## s(d.z):SiteCodePEbrejo  9.00e+00 1.29e+00    1.04    0.99
## s(d.z):SiteCodePEcabo2  9.00e+00 4.49e+00    1.04    0.96
## s(d.z):SiteCodePEcaru1  9.00e+00 5.68e+00    1.04    0.97
## s(d.z):SiteCodePEmata2  9.00e+00 7.83e+00    1.04    0.95
## s(d.z):SiteCodePEsvfer  9.00e+00 8.84e+00    1.04    0.98
## s(d.z):SiteCodePRanto10 9.00e+00 1.00e+00    1.04    0.97
## s(d.z):SiteCodePRanto11 9.00e+00 7.61e+00    1.04    0.97
## s(d.z):SiteCodePRanto12 9.00e+00 3.49e+00    1.04    0.97
## s(d.z):SiteCodePRanto13 9.00e+00 4.82e+00    1.04    0.95
## s(d.z):SiteCodePRanto14 9.00e+00 7.74e+00    1.04    0.99
## s(d.z):SiteCodePRanto15 9.00e+00 3.87e+00    1.04    0.99
## s(d.z):SiteCodePRanto4  9.00e+00 1.00e+00    1.04    0.96
## s(d.z):SiteCodePRanto5  9.00e+00 4.54e+00    1.04    0.98
## s(d.z):SiteCodePRanto6  9.00e+00 3.78e+00    1.04    0.98
## s(d.z):SiteCodePRanto7  9.00e+00 2.43e+00    1.04    0.98
## s(d.z):SiteCodePRanto8  9.00e+00 1.18e+00    1.04    0.98
## s(d.z):SiteCodePRanto9  9.00e+00 1.00e+00    1.04    0.97
## s(d.z):SiteCodePRdiam2  9.00e+00 2.93e+00    1.04    0.97
## s(d.z):SiteCodePRdiam3  9.00e+00 3.51e+00    1.04    0.97
## s(d.z):SiteCodePRdiam4  9.00e+00 6.33e+00    1.04    0.98
## s(d.z):SiteCodePRrico1  9.00e+00 6.55e+00    1.04    0.99
## s(d.z):SiteCodePRrico3  9.00e+00 7.74e+00    1.04    0.98
## s(d.z):SiteCodePRsapo   9.00e+00 8.78e+00    1.04    0.98
## s(d.z):SiteCodePRtele   9.00e+00 7.45e+00    1.04    0.98
## s(d.z):SiteCodePRtiba1  9.00e+00 7.90e+00    1.04    0.97
## s(d.z):SiteCodePRvent   9.00e+00 7.32e+00    1.04    0.98
## s(d.z):SiteCodeRJpnrj   9.00e+00 7.68e+00    1.04    0.95
## s(d.z):SiteCodeRScach1  9.00e+00 6.16e+00    1.04    0.96
## s(d.z):SiteCodeRScach2  9.00e+00 8.10e+00    1.04    0.98
## s(d.z):SiteCodeRScach3  9.00e+00 8.93e+00    1.04    0.97
## s(d.z):SiteCodeRScach4  9.00e+00 4.41e+00    1.04    0.96
## s(d.z):SiteCodeRScama   9.00e+00 7.67e+00    1.04    0.96
## s(d.z):SiteCodeRScamb1  9.00e+00 1.00e+00    1.04    0.97
## s(d.z):SiteCodeRScris   9.00e+00 6.99e+00    1.04    0.97
## s(d.z):SiteCodeRSfaxi   9.00e+00 8.31e+00    1.04    0.98
## s(d.z):SiteCodeRSgaur   9.00e+00 7.26e+00    1.04    0.98
## s(d.z):SiteCodeRSmaqu4  9.00e+00 4.45e+00    1.04    0.98
## s(d.z):SiteCodeRSmorr1  9.00e+00 3.74e+00    1.04    0.98
## s(d.z):SiteCodeRSpet1   9.00e+00 5.52e+00    1.04    0.96
## s(d.z):SiteCodeRSrioz1  9.00e+00 6.84e+00    1.04    0.96
## s(d.z):SiteCodeRSsetape 9.00e+00 7.58e+00    1.04    0.98
## s(d.z):SiteCodeRSsini   9.00e+00 2.56e+00    1.04    0.97
## s(d.z):SiteCodeRSsmar2  9.00e+00 3.14e+00    1.04    0.97
## s(d.z):SiteCodeRSvsol   9.00e+00 7.61e+00    1.04    0.98
## s(d.z):SiteCodeSCarar   9.00e+00 6.02e+00    1.04    0.97
## s(d.z):SiteCodeSCcric1  9.00e+00 7.91e+00    1.04    0.97
## s(d.z):SiteCodeSCilho   9.00e+00 6.28e+00    1.04    0.96
## s(d.z):SiteCodeSCilho2  9.00e+00 4.90e+00    1.04    0.98
## s(d.z):SiteCodeSCorle   9.00e+00 5.29e+00    1.04    0.97
## s(d.z):SiteCodeSCserra4 9.00e+00 1.00e+00    1.04    0.97
## s(d.z):SiteCodeSCside   9.00e+00 1.00e+00    1.04    0.96
## s(d.z):SiteCodeSCside1  9.00e+00 4.10e+00    1.04    0.98
## s(d.z):SiteCodeSCside4  9.00e+00 1.00e+00    1.04    0.99
## s(d.z):SiteCodeSCsjo1   9.00e+00 5.63e+00    1.04    0.98
## s(d.z):SiteCodeSCtimbe1 9.00e+00 2.26e+00    1.04    0.99
## s(d.z):SiteCodeSCvolta1 9.00e+00 7.79e+00    1.04    0.99
## s(d.z):SiteCodeSCvolta2 9.00e+00 4.23e+00    1.04    0.99
## s(d.z):SiteCodeSPcana1  9.00e+00 8.81e+00    1.04    0.96
## s(d.z):SiteCodeSPcara1  9.00e+00 5.15e+00    1.04    0.99
## s(d.z):SiteCodeSPcara3  9.00e+00 8.76e+00    1.04    0.99
## s(d.z):SiteCodeSPeec1   9.00e+00 7.99e+00    1.04    0.99
## s(d.z):SiteCodeSPeec4   9.00e+00 8.34e+00    1.04    0.96
## s(d.z):SiteCodeSPeec5   9.00e+00 7.12e+00    1.04    0.96
## s(d.z):SiteCodeSPigua1  9.00e+00 5.33e+00    1.04    0.98
## s(d.z):SiteCodeSPjabo1  9.00e+00 8.24e+00    1.04    0.97
## s(d.z):SiteCodeSPpecb1  9.00e+00 7.07e+00    1.04    0.97
## s(d.z):SiteCodeSPpei1   9.00e+00 6.97e+00    1.04    0.98
## s(d.z):SiteCodeSPpesm6  9.00e+00 8.54e+00    1.04    0.96
## s(d.z):SiteCodeSPpesmA  9.00e+00 7.99e+00    1.04    0.96
## s(d.z):SiteCodeSPpesmB  9.00e+00 8.05e+00    1.04    0.98
## s(d.z):SiteCodeSPpesmC  9.00e+00 3.27e+00    1.04    0.99
## s(d.z):SiteCodeSPpesmD  9.00e+00 1.00e+00    1.04    0.97
## s(d.z):SiteCodeSPpesmE  9.00e+00 3.89e+00    1.04    0.97
## s(d.z):SiteCodeSPpesmF  9.00e+00 6.67e+00    1.04    0.99
## s(d.z):SiteCodeSPpesmG  9.00e+00 7.21e+00    1.04    0.99
## s(d.z):SiteCodeSPpesmH  9.00e+00 1.83e+00    1.04    0.99
## s(d.z):SiteCodeSPpesmI  9.00e+00 5.21e+00    1.04    0.98
## s(d.z):SiteCodeSPpesmJ  9.00e+00 4.51e+00    1.04    0.95
## s(d.z):SiteCodeSPpesmK  9.00e+00 4.79e+00    1.04    0.96
## s(d.z):SiteCodeSPpesmL  9.00e+00 7.75e+00    1.04    0.99
## s(d.z):SiteCodeSPpesmN  9.00e+00 1.00e+00    1.04    0.98
## s(d.z):SiteCodeSPpesmS  9.00e+00 6.23e+00    1.04    0.99
## s(d.z):SiteCodeSPpesmU  9.00e+00 1.00e+00    1.04    0.99

tabela 2.2.1 Summary

## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## cbind(n_nRef, 100 - n_nRef) ~ s(p.z, bs = "cr") + s(d.z, bs = "cr") + 
##     ti(p.z, d.z, bs = c("cr", "cr")) + s(SiteCode, bs = "re") + 
##     s(d.z, by = SiteCode, bs = "cr", m = 1)
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)    2.985    334.817   0.009    0.993
## 
## Approximate significance of smooth terms:
##                               edf    Ref.df    Chi.sq  p-value    
## s(p.z)                  1.000e+00 1.000e+00     0.000 0.992501    
## s(d.z)                  8.917e+00 8.981e+00   184.902  < 2e-16 ***
## ti(p.z,d.z)             2.217e+00 2.550e+00     5.655 0.090208 .  
## s(SiteCode)             1.007e+02 1.010e+02 16720.810  < 2e-16 ***
## s(d.z):SiteCodeBAjuss   3.518e+00 4.140e+00    55.864 3.48e-11 ***
## s(d.z):SiteCodeBAlenc4  6.213e+00 6.871e+00    22.646 0.001850 ** 
## s(d.z):SiteCodeBAuruc   3.918e+00 4.725e+00    16.256 0.004800 ** 
## s(d.z):SiteCodeESsoor   7.086e+00 7.898e+00    71.275 3.94e-12 ***
## s(d.z):SiteCodeMGipia1  5.600e+00 6.418e+00    19.195 0.005755 ** 
## s(d.z):SiteCodeMGlavr2  7.393e+00 7.876e+00   101.386  < 2e-16 ***
## s(d.z):SiteCodeMGlavr3  1.000e+00 1.000e+00     0.000 0.999850    
## s(d.z):SiteCodeMGlavr6  1.765e+00 1.948e+00     0.122 0.942107    
## s(d.z):SiteCodeMGubera  8.901e+00 8.992e+00   374.208  < 2e-16 ***
## s(d.z):SiteCodeMGuberl1 3.389e+00 3.908e+00    23.179 0.000105 ***
## s(d.z):SiteCodeMGuberl3 2.480e+00 3.080e+00     9.264 0.025249 *  
## s(d.z):SiteCodeMGuberl5 8.906e+00 8.992e+00   128.865  < 2e-16 ***
## s(d.z):SiteCodeMGuberl6 5.254e+00 5.888e+00     9.589 0.115697    
## s(d.z):SiteCodeMGuberl7 4.385e+00 4.918e+00    70.432 1.19e-13 ***
## s(d.z):SiteCodeMGvico1  8.525e+00 8.916e+00   213.180  < 2e-16 ***
## s(d.z):SiteCodeMGvico16 8.150e+00 8.782e+00   118.122  < 2e-16 ***
## s(d.z):SiteCodeMGvico3  7.348e+00 8.071e+00    28.360 0.000512 ***
## s(d.z):SiteCodeMGvico4  8.133e+00 8.698e+00   170.410  < 2e-16 ***
## s(d.z):SiteCodeMSdour   1.628e-06 3.254e-06     0.000 1.000000    
## s(d.z):SiteCodePEalia   6.033e+00 6.761e+00    19.821 0.005221 ** 
## s(d.z):SiteCodePEbrejo  1.294e+00 1.533e+00     0.263 0.782471    
## s(d.z):SiteCodePEcabo2  4.492e+00 5.453e+00    10.946 0.064084 .  
## s(d.z):SiteCodePEcaru1  5.676e+00 6.434e+00     8.413 0.236611    
## s(d.z):SiteCodePEmata2  7.827e+00 8.582e+00   114.401  < 2e-16 ***
## s(d.z):SiteCodePEsvfer  8.837e+00 8.974e+00     1.257 0.998538    
## s(d.z):SiteCodePRanto10 1.000e+00 1.001e+00     0.121 0.728520    
## s(d.z):SiteCodePRanto11 7.611e+00 8.318e+00    32.710 8.71e-05 ***
## s(d.z):SiteCodePRanto12 3.490e+00 4.185e+00     6.615 0.170366    
## s(d.z):SiteCodePRanto13 4.820e+00 5.781e+00    21.908 0.001183 ** 
## s(d.z):SiteCodePRanto14 7.735e+00 8.369e+00    36.584 1.58e-05 ***
## s(d.z):SiteCodePRanto15 3.869e+00 4.386e+00     6.884 0.215468    
## s(d.z):SiteCodePRanto4  1.000e+00 1.000e+00     0.118 0.731777    
## s(d.z):SiteCodePRanto5  4.537e+00 5.491e+00    15.817 0.010771 *  
## s(d.z):SiteCodePRanto6  3.780e+00 4.453e+00    30.756 7.98e-06 ***
## s(d.z):SiteCodePRanto7  2.431e+00 2.915e+00    13.806 0.006562 ** 
## s(d.z):SiteCodePRanto8  1.181e+00 1.789e+00     1.340 0.496320    
## s(d.z):SiteCodePRanto9  1.000e+00 1.000e+00     0.130 0.718262    
## s(d.z):SiteCodePRdiam2  2.929e+00 3.334e+00     7.946 0.070747 .  
## s(d.z):SiteCodePRdiam3  3.509e+00 4.042e+00    37.371 1.60e-07 ***
## s(d.z):SiteCodePRdiam4  6.334e+00 7.126e+00     8.110 0.334557    
## s(d.z):SiteCodePRrico1  6.553e+00 7.391e+00    98.485  < 2e-16 ***
## s(d.z):SiteCodePRrico3  7.743e+00 8.230e+00    81.085 2.30e-13 ***
## s(d.z):SiteCodePRsapo   8.782e+00 8.967e+00    94.450  < 2e-16 ***
## s(d.z):SiteCodePRtele   7.446e+00 8.185e+00   127.443  < 2e-16 ***
## s(d.z):SiteCodePRtiba1  7.895e+00 8.546e+00   130.425  < 2e-16 ***
## s(d.z):SiteCodePRvent   7.316e+00 7.892e+00   137.352  < 2e-16 ***
## s(d.z):SiteCodeRJpnrj   7.676e+00 8.333e+00    71.624 4.39e-12 ***
## s(d.z):SiteCodeRScach1  6.162e+00 6.781e+00    52.604 7.60e-09 ***
## s(d.z):SiteCodeRScach2  8.105e+00 8.620e+00   156.415  < 2e-16 ***
## s(d.z):SiteCodeRScach3  8.930e+00 8.997e+00   214.992  < 2e-16 ***
## s(d.z):SiteCodeRScach4  4.410e+00 4.998e+00     9.066 0.113431    
## s(d.z):SiteCodeRScama   7.671e+00 8.277e+00    17.616 0.028409 *  
## s(d.z):SiteCodeRScamb1  1.000e+00 1.000e+00     0.000 0.995301    
## s(d.z):SiteCodeRScris   6.991e+00 7.624e+00    45.902 6.73e-07 ***
## s(d.z):SiteCodeRSfaxi   8.312e+00 8.760e+00   172.834  < 2e-16 ***
## s(d.z):SiteCodeRSgaur   7.264e+00 7.919e+00    31.772 9.40e-05 ***
## s(d.z):SiteCodeRSmaqu4  4.452e+00 5.426e+00    44.325 8.46e-08 ***
## s(d.z):SiteCodeRSmorr1  3.741e+00 4.418e+00     7.862 0.111616    
## s(d.z):SiteCodeRSpet1   5.521e+00 6.244e+00    50.795 6.03e-09 ***
## s(d.z):SiteCodeRSrioz1  6.841e+00 7.188e+00     0.515 0.999512    
## s(d.z):SiteCodeRSsetape 7.577e+00 7.892e+00    14.478 0.074150 .  
## s(d.z):SiteCodeRSsini   2.559e+00 3.215e+00     5.872 0.091947 .  
## s(d.z):SiteCodeRSsmar2  3.143e+00 3.599e+00    12.399 0.019760 *  
## s(d.z):SiteCodeRSvsol   7.609e+00 8.296e+00    43.975 7.49e-07 ***
## s(d.z):SiteCodeSCarar   6.021e+00 6.762e+00     7.108 0.390389    
## s(d.z):SiteCodeSCcric1  7.914e+00 8.496e+00   116.890  < 2e-16 ***
## s(d.z):SiteCodeSCilho   6.284e+00 7.056e+00    49.581 1.92e-08 ***
## s(d.z):SiteCodeSCilho2  4.900e+00 5.863e+00    17.084 0.008020 ** 
## s(d.z):SiteCodeSCorle   5.294e+00 5.952e+00    77.425 1.17e-14 ***
## s(d.z):SiteCodeSCserra4 1.000e+00 1.001e+00     0.000 0.993642    
## s(d.z):SiteCodeSCside   1.000e+00 1.001e+00     0.000 0.990357    
## s(d.z):SiteCodeSCside1  4.103e+00 4.842e+00    38.398 3.26e-07 ***
## s(d.z):SiteCodeSCside4  1.000e+00 1.000e+00     0.000 0.996727    
## s(d.z):SiteCodeSCsjo1   5.627e+00 6.282e+00    27.456 0.000167 ***
## s(d.z):SiteCodeSCtimbe1 2.263e+00 2.833e+00     5.977 0.177634    
## s(d.z):SiteCodeSCvolta1 7.792e+00 8.490e+00   102.101  < 2e-16 ***
## s(d.z):SiteCodeSCvolta2 4.227e+00 5.001e+00    18.193 0.002718 ** 
## s(d.z):SiteCodeSPcana1  8.810e+00 8.981e+00   108.511  < 2e-16 ***
## s(d.z):SiteCodeSPcara1  5.154e+00 6.079e+00    17.558 0.007704 ** 
## s(d.z):SiteCodeSPcara3  8.763e+00 8.966e+00   126.650  < 2e-16 ***
## s(d.z):SiteCodeSPeec1   7.991e+00 8.609e+00   132.468  < 2e-16 ***
## s(d.z):SiteCodeSPeec4   8.341e+00 8.786e+00    25.417 0.001971 ** 
## s(d.z):SiteCodeSPeec5   7.123e+00 7.828e+00    45.583 3.12e-07 ***
## s(d.z):SiteCodeSPigua1  5.332e+00 6.161e+00    14.986 0.025390 *  
## s(d.z):SiteCodeSPjabo1  8.242e+00 8.718e+00   255.351  < 2e-16 ***
## s(d.z):SiteCodeSPpecb1  7.066e+00 7.735e+00    50.802 2.18e-08 ***
## s(d.z):SiteCodeSPpei1   6.966e+00 7.873e+00    74.614 6.27e-13 ***
## s(d.z):SiteCodeSPpesm6  8.537e+00 8.898e+00   114.199  < 2e-16 ***
## s(d.z):SiteCodeSPpesmA  7.988e+00 8.645e+00    41.836 3.61e-06 ***
## s(d.z):SiteCodeSPpesmB  8.046e+00 8.598e+00    14.627 0.139066    
## s(d.z):SiteCodeSPpesmC  3.271e+00 3.928e+00     4.808 0.239948    
## s(d.z):SiteCodeSPpesmD  1.000e+00 1.001e+00     0.000 0.993652    
## s(d.z):SiteCodeSPpesmE  3.895e+00 4.591e+00     6.616 0.182925    
## s(d.z):SiteCodeSPpesmF  6.674e+00 7.643e+00    11.998 0.128377    
## s(d.z):SiteCodeSPpesmG  7.208e+00 7.967e+00    20.811 0.009023 ** 
## s(d.z):SiteCodeSPpesmH  1.825e+00 2.285e+00     1.549 0.460585    
## s(d.z):SiteCodeSPpesmI  5.206e+00 6.062e+00    20.324 0.002539 ** 
## s(d.z):SiteCodeSPpesmJ  4.505e+00 5.275e+00     8.113 0.173980    
## s(d.z):SiteCodeSPpesmK  4.786e+00 5.615e+00    27.558 6.19e-05 ***
## s(d.z):SiteCodeSPpesmL  7.750e+00 8.423e+00    38.378 9.06e-06 ***
## s(d.z):SiteCodeSPpesmN  1.000e+00 1.001e+00     0.007 0.932461    
## s(d.z):SiteCodeSPpesmS  6.230e+00 7.199e+00     8.561 0.292272    
## s(d.z):SiteCodeSPpesmU  1.000e+00 1.000e+00     0.000 0.985334    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 1063/1065
## R-sq.(adj) =  0.986   Deviance explained = 98.3%
## UBRE = 0.69728  Scale est. = 1         n = 2060

Figura 2.2.2 Observado e predito pelo modelo cheio mais plausível

Figura 2.2.3 Diagnostico: avaliação do ajuste pelas preditoras lineares

v_sitesOut_nRefEE <- filter(df_plot,.fitted < -300 | .fitted > 380 | .cooksd > 2000) %>%
  .$SiteCode %>% unique

Comparação Modelo Cheio sem outliers

# sites outlier
v_sitesOut_nRefEE <- filter(df_plot,.fitted < -300 | .fitted > 380 | .cooksd > 2000) %>%
  .$SiteCode %>% unique
# dados
df_md <- df_resultados %>% filter(MN=="EE" & !(SiteCode %in% v_sitesOut_nRefEE)) %>% 
  distinct() %>% 
  mutate(rep = 1) 
df_md <- left_join(x=data.frame(rep=1,
                                md_class=c("glmm d+d^2|Site",
                                           "glmm d|Site",
                                           "glmm 1|Site",
                                           "gamm d|Site for each",
                                           "gamm d|Site common",
                                           "gamm 1|Site")),
                   y=df_md,
                   by="rep")
# ajuste modelos
registerDoMC(3)
l_nRefEE_mdCheio__sOut <- dlply(df_md,"md_class",f_md,.parallel = TRUE)
##                      dAICc   df               weight
## gamm d|Site for each     0.0 644.383383100452 1     
## gamm d|Site common     778.8 826.614698460209 <0.001
## glmm d+d^2|Site       5111.4 15               <0.001
## glmm d|Site           8063.3 12               <0.001
## gamm 1|Site          23983.9 123.849903818296 <0.001
## glmm 1|Site          26005.3 10               <0.001

Diagnosticos modelo cheio sem outliers

Figura 2.2.4 Appraise(md_nRefEE_sOut)

Figura 2.2.5 Diagnostico: avaliação do ajuste pelas preditoras lineares

Subset=MNEI

Estudo GAMM: smoother type

  1. Antes de comparar os modelos cheios preciso avaliar smoother type mais para utilizar na GAMM;
  2. Pelo observado nos dados (figuras 2.1 e 2.2) vou comparar dois smoother type: “tp” e “cr”
  3. Para isso vou ajuster um smoother para d por Sítio de amostragem com um mesmo parâmetro de penalização, dessa maneira eu considero a variação por sítio de amostragem na variável d sem um elevado custo computacional.
# dados
df_md <- df_resultados %>% filter(MN=="EI") %>% distinct()
df_md <- rbind(mutate(df_md,model_class = "cr"),
               mutate(df_md,model_class = "tp"))
df_md$model_class <- factor(df_md$model_class)
# função
f_md <- function(dados){
  model_class <- unique(dados$model_class)
  if(model_class == "cr"){
    md_ <- gam(cbind(n_nRef,100-n_nRef) ~ 
                 s(p.z,bs="cr") + s(d.z,bs="cr") + ti(p.z,d.z,bs=c("cr","cr")) +
                 s(d.z,SiteCode,bs="fs",xt=list(bs="cr"), m=1),
                data = dados, family = "binomial")
    
  }else{
    md_ <- gam(cbind(n_nRef,100-n_nRef) ~ 
                 s(p.z,bs="tp") + s(d.z,bs="tp") + ti(p.z,d.z,bs=c("tp","tp")) +
                 s(d.z,SiteCode,bs="fs",xt=list(bs="tp"), m=1),
                data = dados, family = "binomial")
    
  }
  return(md_)
}
registerDoMC(2)
l_nRefEI_estudoGAMM <- dlply(df_md,"model_class",f_md,.parallel = TRUE)
##    dAICc df               weight
## cr   0.0 647.199595731505 1     
## tp 443.6 758.152972168071 <0.001

Comparação Modelos Cheios

# dados
df_md <- df_resultados %>% filter(MN=="EI") %>% 
  distinct() %>% 
  mutate(rep = 1)
df_md <- left_join(x=data.frame(rep=1,
                                md_class=c("glmm d+d^2|Site",
                                           "glmm d|Site",
                                           "glmm 1|Site",
                                           "gamm d|Site for each",
                                           "gamm d|Site common",
                                           "gamm 1|Site")),
                   y=df_md,
                   by="rep")
# função de ajuste
registerDoMC(3)
l_nRefEI_mdCheio <- dlply(df_md,"md_class",f_md,.parallel = TRUE)
##                      dAICc   df               weight
## gamm d|Site for each     0.0 556.021191967948 1     
## gamm d|Site common     374.1 647.199595731505 <0.001
## glmm d+d^2|Site       3256.4 15               <0.001
## glmm d|Site          11864.4 12               <0.001
## gamm 1|Site          57181.1 125.994597023309 <0.001
## glmm 1|Site          63995.3 10               <0.001

Diagnostico modelo cheio

Figura 2.3.1 appraise(md_nRefEI)

## 
## Method: UBRE   Optimizer: outer newton
## full convergence after 50 iterations.
## Gradient range [-2.37292e-06,8.788319e-07]
## (score 0.1963635 & scale 1).
## eigenvalue range [-7.886026e-07,0.001085589].
## Model rank =  1063 / 1065 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                               k'      edf k-index p-value   
## s(p.z)                  9.00e+00 1.00e+00    1.25   1.000   
## s(d.z)                  9.00e+00 7.25e+00    0.96   0.020 * 
## ti(p.z,d.z)             1.60e+01 1.29e+01    1.00   0.640   
## s(SiteCode)             1.03e+02 1.01e+02      NA      NA   
## s(d.z):SiteCodeBAjuss   9.00e+00 2.88e+00    0.96   0.040 * 
## s(d.z):SiteCodeBAlenc4  9.00e+00 3.45e+00    0.96   0.025 * 
## s(d.z):SiteCodeBAuruc   9.00e+00 1.00e+00    0.96   0.045 * 
## s(d.z):SiteCodeESsoor   9.00e+00 2.52e+00    0.96   0.015 * 
## s(d.z):SiteCodeMGipia1  9.00e+00 3.53e+00    0.96   0.065 . 
## s(d.z):SiteCodeMGlavr2  9.00e+00 8.56e+00    0.96   0.035 * 
## s(d.z):SiteCodeMGlavr3  9.00e+00 5.20e-05    0.96   0.020 * 
## s(d.z):SiteCodeMGlavr6  9.00e+00 7.22e+00    0.96   0.015 * 
## s(d.z):SiteCodeMGubera  9.00e+00 6.97e+00    0.96   0.040 * 
## s(d.z):SiteCodeMGuberl1 9.00e+00 4.87e+00    0.96   0.050 * 
## s(d.z):SiteCodeMGuberl3 9.00e+00 4.06e+00    0.96   0.015 * 
## s(d.z):SiteCodeMGuberl5 9.00e+00 5.98e+00    0.96   0.020 * 
## s(d.z):SiteCodeMGuberl6 9.00e+00 4.05e+00    0.96   0.025 * 
## s(d.z):SiteCodeMGuberl7 9.00e+00 3.50e+00    0.96   0.045 * 
## s(d.z):SiteCodeMGvico1  9.00e+00 1.00e+00    0.96   0.025 * 
## s(d.z):SiteCodeMGvico16 9.00e+00 1.00e+00    0.96   0.045 * 
## s(d.z):SiteCodeMGvico3  9.00e+00 1.07e+00    0.96   0.050 * 
## s(d.z):SiteCodeMGvico4  9.00e+00 1.00e+00    0.96   0.025 * 
## s(d.z):SiteCodeMSdour   9.00e+00 2.53e+00    0.96   0.035 * 
## s(d.z):SiteCodePEalia   9.00e+00 2.41e+00    0.96   0.045 * 
## s(d.z):SiteCodePEbrejo  9.00e+00 5.62e+00    0.96   0.015 * 
## s(d.z):SiteCodePEcabo2  9.00e+00 4.30e+00    0.96   0.045 * 
## s(d.z):SiteCodePEcaru1  9.00e+00 5.14e+00    0.96   0.070 . 
## s(d.z):SiteCodePEmata2  9.00e+00 1.00e+00    0.96   0.045 * 
## s(d.z):SiteCodePEsvfer  9.00e+00 8.26e+00    0.96   0.050 * 
## s(d.z):SiteCodePRanto10 9.00e+00 6.33e+00    0.96   0.045 * 
## s(d.z):SiteCodePRanto11 9.00e+00 3.68e+00    0.96   0.055 . 
## s(d.z):SiteCodePRanto12 9.00e+00 7.30e+00    0.96   0.040 * 
## s(d.z):SiteCodePRanto13 9.00e+00 1.00e+00    0.96   0.025 * 
## s(d.z):SiteCodePRanto14 9.00e+00 5.90e+00    0.96   0.040 * 
## s(d.z):SiteCodePRanto15 9.00e+00 2.85e+00    0.96   0.040 * 
## s(d.z):SiteCodePRanto4  9.00e+00 5.53e+00    0.96   0.040 * 
## s(d.z):SiteCodePRanto5  9.00e+00 1.00e+00    0.96   0.030 * 
## s(d.z):SiteCodePRanto6  9.00e+00 6.85e+00    0.96   0.045 * 
## s(d.z):SiteCodePRanto7  9.00e+00 5.64e+00    0.96   0.050 * 
## s(d.z):SiteCodePRanto8  9.00e+00 1.00e+00    0.96   0.015 * 
## s(d.z):SiteCodePRanto9  9.00e+00 5.01e+00    0.96   0.035 * 
## s(d.z):SiteCodePRdiam2  9.00e+00 1.23e+00    0.96   0.020 * 
## s(d.z):SiteCodePRdiam3  9.00e+00 1.02e+00    0.96   0.040 * 
## s(d.z):SiteCodePRdiam4  9.00e+00 2.88e+00    0.96   0.020 * 
## s(d.z):SiteCodePRrico1  9.00e+00 4.73e+00    0.96   0.025 * 
## s(d.z):SiteCodePRrico3  9.00e+00 1.00e+00    0.96   0.045 * 
## s(d.z):SiteCodePRsapo   9.00e+00 6.54e+00    0.96   0.050 * 
## s(d.z):SiteCodePRtele   9.00e+00 5.02e+00    0.96   0.030 * 
## s(d.z):SiteCodePRtiba1  9.00e+00 2.65e+00    0.96   0.045 * 
## s(d.z):SiteCodePRvent   9.00e+00 5.25e+00    0.96   0.035 * 
## s(d.z):SiteCodeRJpnrj   9.00e+00 1.00e+00    0.96   0.070 . 
## s(d.z):SiteCodeRScach1  9.00e+00 3.56e+00    0.96   0.040 * 
## s(d.z):SiteCodeRScach2  9.00e+00 2.39e+00    0.96   0.045 * 
## s(d.z):SiteCodeRScach3  9.00e+00 2.83e+00    0.96   0.015 * 
## s(d.z):SiteCodeRScach4  9.00e+00 1.00e+00    0.96   0.060 . 
## s(d.z):SiteCodeRScama   9.00e+00 7.58e+00    0.96   0.015 * 
## s(d.z):SiteCodeRScamb1  9.00e+00 3.01e+00    0.96   0.030 * 
## s(d.z):SiteCodeRScris   9.00e+00 5.40e+00    0.96   0.065 . 
## s(d.z):SiteCodeRSfaxi   9.00e+00 4.59e+00    0.96   0.035 * 
## s(d.z):SiteCodeRSgaur   9.00e+00 4.13e+00    0.96   0.020 * 
## s(d.z):SiteCodeRSmaqu4  9.00e+00 7.67e+00    0.96   0.040 * 
## s(d.z):SiteCodeRSmorr1  9.00e+00 6.91e+00    0.96   0.035 * 
## s(d.z):SiteCodeRSpet1   9.00e+00 5.21e+00    0.96   0.055 . 
## s(d.z):SiteCodeRSrioz1  9.00e+00 5.08e+00    0.96   0.040 * 
## s(d.z):SiteCodeRSsetape 9.00e+00 5.86e+00    0.96   0.030 * 
## s(d.z):SiteCodeRSsini   9.00e+00 6.37e+00    0.96   0.035 * 
## s(d.z):SiteCodeRSsmar2  9.00e+00 8.17e+00    0.96   0.055 . 
## s(d.z):SiteCodeRSvsol   9.00e+00 2.73e+00    0.96   0.010 **
## s(d.z):SiteCodeSCarar   9.00e+00 1.51e+00    0.96   0.030 * 
## s(d.z):SiteCodeSCcric1  9.00e+00 3.16e+00    0.96   0.025 * 
## s(d.z):SiteCodeSCilho   9.00e+00 2.98e+00    0.96   0.045 * 
## s(d.z):SiteCodeSCilho2  9.00e+00 4.10e+00    0.96   0.055 . 
## s(d.z):SiteCodeSCorle   9.00e+00 5.77e+00    0.96   0.040 * 
## s(d.z):SiteCodeSCserra4 9.00e+00 6.62e+00    0.96   0.045 * 
## s(d.z):SiteCodeSCside   9.00e+00 3.28e+00    0.96   0.045 * 
## s(d.z):SiteCodeSCside1  9.00e+00 2.71e+00    0.96   0.025 * 
## s(d.z):SiteCodeSCside4  9.00e+00 1.00e+00    0.96   0.020 * 
## s(d.z):SiteCodeSCsjo1   9.00e+00 4.88e+00    0.96   0.060 . 
## s(d.z):SiteCodeSCtimbe1 9.00e+00 6.09e+00    0.96   0.035 * 
## s(d.z):SiteCodeSCvolta1 9.00e+00 4.50e+00    0.96   0.030 * 
## s(d.z):SiteCodeSCvolta2 9.00e+00 5.90e+00    0.96   0.020 * 
## s(d.z):SiteCodeSPcana1  9.00e+00 5.20e+00    0.96   0.030 * 
## s(d.z):SiteCodeSPcara1  9.00e+00 7.14e+00    0.96   0.030 * 
## s(d.z):SiteCodeSPcara3  9.00e+00 7.60e+00    0.96   0.045 * 
## s(d.z):SiteCodeSPeec1   9.00e+00 1.00e+00    0.96   0.055 . 
## s(d.z):SiteCodeSPeec4   9.00e+00 6.21e+00    0.96   0.035 * 
## s(d.z):SiteCodeSPeec5   9.00e+00 1.59e+00    0.96   0.045 * 
## s(d.z):SiteCodeSPigua1  9.00e+00 1.00e+00    0.96   0.055 . 
## s(d.z):SiteCodeSPjabo1  9.00e+00 6.59e+00    0.96   0.025 * 
## s(d.z):SiteCodeSPpecb1  9.00e+00 1.93e+00    0.96   0.020 * 
## s(d.z):SiteCodeSPpei1   9.00e+00 5.09e+00    0.96   0.050 * 
## s(d.z):SiteCodeSPpesm6  9.00e+00 3.56e+00    0.96   0.035 * 
## s(d.z):SiteCodeSPpesmA  9.00e+00 5.78e+00    0.96   0.025 * 
## s(d.z):SiteCodeSPpesmB  9.00e+00 3.44e+00    0.96   0.015 * 
## s(d.z):SiteCodeSPpesmC  9.00e+00 3.10e+00    0.96   0.020 * 
## s(d.z):SiteCodeSPpesmD  9.00e+00 3.66e+00    0.96   0.030 * 
## s(d.z):SiteCodeSPpesmE  9.00e+00 4.04e+00    0.96   0.035 * 
## s(d.z):SiteCodeSPpesmF  9.00e+00 3.50e+00    0.96   0.065 . 
## s(d.z):SiteCodeSPpesmG  9.00e+00 5.09e+00    0.96   0.045 * 
## s(d.z):SiteCodeSPpesmH  9.00e+00 5.30e+00    0.96   0.045 * 
## s(d.z):SiteCodeSPpesmI  9.00e+00 6.60e+00    0.96   0.035 * 
## s(d.z):SiteCodeSPpesmJ  9.00e+00 6.83e+00    0.96   0.045 * 
## s(d.z):SiteCodeSPpesmK  9.00e+00 4.66e+00    0.96   0.050 * 
## s(d.z):SiteCodeSPpesmL  9.00e+00 5.60e+00    0.96   0.025 * 
## s(d.z):SiteCodeSPpesmN  9.00e+00 4.81e+00    0.96   0.025 * 
## s(d.z):SiteCodeSPpesmS  9.00e+00 5.19e+00    0.96   0.025 * 
## s(d.z):SiteCodeSPpesmU  9.00e+00 8.32e+00    0.96   0.025 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## cbind(n_nRef, 100 - n_nRef) ~ s(p.z, bs = "cr") + s(d.z, bs = "cr") + 
##     ti(p.z, d.z, bs = c("cr", "cr")) + s(SiteCode, bs = "re") + 
##     s(d.z, by = SiteCode, bs = "cr", m = 1)
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  -0.5485    10.6261  -0.052    0.959
## 
## Approximate significance of smooth terms:
##                               edf    Ref.df    Chi.sq  p-value    
## s(p.z)                  1.000e+00 1.000e+00     0.001 0.969635    
## s(d.z)                  7.249e+00 7.995e+00   599.870  < 2e-16 ***
## ti(p.z,d.z)             1.288e+01 1.296e+01  1918.523  < 2e-16 ***
## s(SiteCode)             1.008e+02 1.010e+02 12589.462  < 2e-16 ***
## s(d.z):SiteCodeBAjuss   2.885e+00 3.083e+00     7.211 0.043988 *  
## s(d.z):SiteCodeBAlenc4  3.450e+00 4.178e+00    14.025 0.008454 ** 
## s(d.z):SiteCodeBAuruc   1.000e+00 1.000e+00     2.169 0.140871    
## s(d.z):SiteCodeESsoor   2.520e+00 2.947e+00    11.587 0.016453 *  
## s(d.z):SiteCodeMGipia1  3.526e+00 4.180e+00    10.011 0.042309 *  
## s(d.z):SiteCodeMGlavr2  8.565e+00 8.914e+00    33.426 0.000297 ***
## s(d.z):SiteCodeMGlavr3  5.202e-05 9.743e-05     0.000 0.999819    
## s(d.z):SiteCodeMGlavr6  7.224e+00 8.001e+00    19.314 0.013379 *  
## s(d.z):SiteCodeMGubera  6.967e+00 7.543e+00    88.994 1.27e-15 ***
## s(d.z):SiteCodeMGuberl1 4.869e+00 5.678e+00   152.128  < 2e-16 ***
## s(d.z):SiteCodeMGuberl3 4.057e+00 4.755e+00    50.761 1.09e-09 ***
## s(d.z):SiteCodeMGuberl5 5.977e+00 6.744e+00    41.273 7.68e-07 ***
## s(d.z):SiteCodeMGuberl6 4.053e+00 4.746e+00    69.693 1.86e-13 ***
## s(d.z):SiteCodeMGuberl7 3.503e+00 4.176e+00    22.220 0.000249 ***
## s(d.z):SiteCodeMGvico1  1.000e+00 1.000e+00    10.331 0.001310 ** 
## s(d.z):SiteCodeMGvico16 1.000e+00 1.000e+00    10.693 0.001075 ** 
## s(d.z):SiteCodeMGvico3  1.066e+00 1.127e+00    11.695 0.001110 ** 
## s(d.z):SiteCodeMGvico4  1.002e+00 1.004e+00     9.709 0.001845 ** 
## s(d.z):SiteCodeMSdour   2.529e+00 3.209e+00    16.050 0.001154 ** 
## s(d.z):SiteCodePEalia   2.409e+00 2.912e+00    24.967 0.000137 ***
## s(d.z):SiteCodePEbrejo  5.618e+00 6.529e+00   102.805  < 2e-16 ***
## s(d.z):SiteCodePEcabo2  4.296e+00 5.148e+00    96.974  < 2e-16 ***
## s(d.z):SiteCodePEcaru1  5.137e+00 5.974e+00   125.181  < 2e-16 ***
## s(d.z):SiteCodePEmata2  1.000e+00 1.000e+00     2.394 0.121837    
## s(d.z):SiteCodePEsvfer  8.264e+00 8.806e+00   162.156  < 2e-16 ***
## s(d.z):SiteCodePRanto10 6.330e+00 7.239e+00   126.457  < 2e-16 ***
## s(d.z):SiteCodePRanto11 3.679e+00 4.402e+00    27.440 6.02e-05 ***
## s(d.z):SiteCodePRanto12 7.299e+00 8.174e+00    60.495 4.90e-10 ***
## s(d.z):SiteCodePRanto13 1.001e+00 1.002e+00     1.105 0.294043    
## s(d.z):SiteCodePRanto14 5.902e+00 7.076e+00    37.792 3.99e-06 ***
## s(d.z):SiteCodePRanto15 2.853e+00 3.519e+00     5.771 0.315062    
## s(d.z):SiteCodePRanto4  5.527e+00 6.364e+00    66.910 3.70e-12 ***
## s(d.z):SiteCodePRanto5  1.000e+00 1.001e+00     1.044 0.307079    
## s(d.z):SiteCodePRanto6  6.845e+00 7.812e+00     9.459 0.291539    
## s(d.z):SiteCodePRanto7  5.640e+00 6.689e+00   359.758  < 2e-16 ***
## s(d.z):SiteCodePRanto8  1.004e+00 1.009e+00     1.117 0.291522    
## s(d.z):SiteCodePRanto9  5.009e+00 5.897e+00   144.199  < 2e-16 ***
## s(d.z):SiteCodePRdiam2  1.226e+00 1.418e+00    10.480 0.005616 ** 
## s(d.z):SiteCodePRdiam3  1.015e+00 1.030e+00     8.877 0.002987 ** 
## s(d.z):SiteCodePRdiam4  2.878e+00 3.612e+00    13.345 0.010096 *  
## s(d.z):SiteCodePRrico1  4.728e+00 5.406e+00    14.379 0.019964 *  
## s(d.z):SiteCodePRrico3  1.000e+00 1.000e+00     7.989 0.004707 ** 
## s(d.z):SiteCodePRsapo   6.541e+00 7.501e+00   330.201  < 2e-16 ***
## s(d.z):SiteCodePRtele   5.021e+00 5.832e+00    17.533 0.005833 ** 
## s(d.z):SiteCodePRtiba1  2.646e+00 3.260e+00    10.179 0.019416 *  
## s(d.z):SiteCodePRvent   5.248e+00 6.146e+00    94.939  < 2e-16 ***
## s(d.z):SiteCodeRJpnrj   1.000e+00 1.000e+00     7.511 0.006136 ** 
## s(d.z):SiteCodeRScach1  3.557e+00 4.338e+00    34.071 1.13e-06 ***
## s(d.z):SiteCodeRScach2  2.389e+00 3.055e+00    12.198 0.006591 ** 
## s(d.z):SiteCodeRScach3  2.833e+00 3.439e+00    63.362 6.16e-13 ***
## s(d.z):SiteCodeRScach4  1.000e+00 1.000e+00     9.436 0.002128 ** 
## s(d.z):SiteCodeRScama   7.584e+00 8.323e+00   209.383  < 2e-16 ***
## s(d.z):SiteCodeRScamb1  3.008e+00 3.665e+00    52.535 7.17e-11 ***
## s(d.z):SiteCodeRScris   5.399e+00 6.118e+00    58.011 1.41e-10 ***
## s(d.z):SiteCodeRSfaxi   4.590e+00 5.314e+00   307.812  < 2e-16 ***
## s(d.z):SiteCodeRSgaur   4.128e+00 4.885e+00   286.344  < 2e-16 ***
## s(d.z):SiteCodeRSmaqu4  7.666e+00 8.415e+00    45.417 4.66e-07 ***
## s(d.z):SiteCodeRSmorr1  6.906e+00 7.821e+00   293.708  < 2e-16 ***
## s(d.z):SiteCodeRSpet1   5.209e+00 6.231e+00    10.955 0.074749 .  
## s(d.z):SiteCodeRSrioz1  5.080e+00 5.964e+00    31.943 1.92e-05 ***
## s(d.z):SiteCodeRSsetape 5.856e+00 6.568e+00   167.875  < 2e-16 ***
## s(d.z):SiteCodeRSsini   6.372e+00 7.146e+00    72.276 6.73e-13 ***
## s(d.z):SiteCodeRSsmar2  8.174e+00 8.716e+00   244.028  < 2e-16 ***
## s(d.z):SiteCodeRSvsol   2.727e+00 3.283e+00    51.506 9.54e-10 ***
## s(d.z):SiteCodeSCarar   1.505e+00 1.871e+00    13.075 0.003705 ** 
## s(d.z):SiteCodeSCcric1  3.163e+00 3.925e+00    17.572 0.001897 ** 
## s(d.z):SiteCodeSCilho   2.976e+00 3.672e+00    23.916 0.000124 ***
## s(d.z):SiteCodeSCilho2  4.097e+00 4.893e+00    51.290 6.65e-10 ***
## s(d.z):SiteCodeSCorle   5.769e+00 6.826e+00   211.219  < 2e-16 ***
## s(d.z):SiteCodeSCserra4 6.620e+00 7.540e+00    81.477 1.59e-14 ***
## s(d.z):SiteCodeSCside   3.278e+00 4.059e+00    16.239 0.002868 ** 
## s(d.z):SiteCodeSCside1  2.710e+00 3.297e+00    76.515 5.85e-16 ***
## s(d.z):SiteCodeSCside4  1.000e+00 1.000e+00     0.072 0.788352    
## s(d.z):SiteCodeSCsjo1   4.882e+00 5.747e+00   226.634  < 2e-16 ***
## s(d.z):SiteCodeSCtimbe1 6.090e+00 7.086e+00   301.212  < 2e-16 ***
## s(d.z):SiteCodeSCvolta1 4.497e+00 5.367e+00   144.630  < 2e-16 ***
## s(d.z):SiteCodeSCvolta2 5.898e+00 6.835e+00   164.919  < 2e-16 ***
## s(d.z):SiteCodeSPcana1  5.196e+00 6.341e+00    19.944 0.003461 ** 
## s(d.z):SiteCodeSPcara1  7.135e+00 8.026e+00   326.287  < 2e-16 ***
## s(d.z):SiteCodeSPcara3  7.604e+00 8.432e+00   312.611  < 2e-16 ***
## s(d.z):SiteCodeSPeec1   1.001e+00 1.002e+00    11.864 0.000575 ***
## s(d.z):SiteCodeSPeec4   6.214e+00 7.056e+00    15.967 0.026046 *  
## s(d.z):SiteCodeSPeec5   1.585e+00 1.936e+00     9.642 0.005023 ** 
## s(d.z):SiteCodeSPigua1  1.000e+00 1.000e+00     3.236 0.072067 .  
## s(d.z):SiteCodeSPjabo1  6.590e+00 7.380e+00    20.934 0.006802 ** 
## s(d.z):SiteCodeSPpecb1  1.930e+00 2.544e+00    19.084 0.002274 ** 
## s(d.z):SiteCodeSPpei1   5.087e+00 6.116e+00     9.397 0.168178    
## s(d.z):SiteCodeSPpesm6  3.563e+00 4.371e+00   155.077  < 2e-16 ***
## s(d.z):SiteCodeSPpesmA  5.780e+00 6.802e+00   129.630  < 2e-16 ***
## s(d.z):SiteCodeSPpesmB  3.443e+00 4.133e+00    89.219  < 2e-16 ***
## s(d.z):SiteCodeSPpesmC  3.098e+00 3.763e+00    26.642 2.77e-05 ***
## s(d.z):SiteCodeSPpesmD  3.661e+00 4.266e+00   117.521  < 2e-16 ***
## s(d.z):SiteCodeSPpesmE  4.039e+00 4.853e+00    39.450 2.11e-07 ***
## s(d.z):SiteCodeSPpesmF  3.499e+00 4.188e+00    36.301 3.88e-07 ***
## s(d.z):SiteCodeSPpesmG  5.087e+00 6.079e+00    14.182 0.027445 *  
## s(d.z):SiteCodeSPpesmH  5.305e+00 6.211e+00    16.330 0.013242 *  
## s(d.z):SiteCodeSPpesmI  6.596e+00 7.549e+00    64.524 7.11e-11 ***
## s(d.z):SiteCodeSPpesmJ  6.829e+00 7.609e+00    66.533 2.29e-11 ***
## s(d.z):SiteCodeSPpesmK  4.660e+00 5.463e+00    53.814 7.52e-10 ***
## s(d.z):SiteCodeSPpesmL  5.601e+00 6.492e+00    86.200 7.25e-16 ***
## s(d.z):SiteCodeSPpesmN  4.807e+00 5.718e+00    37.147 1.12e-06 ***
## s(d.z):SiteCodeSPpesmS  5.188e+00 6.096e+00   281.815  < 2e-16 ***
## s(d.z):SiteCodeSPpesmU  8.322e+00 8.808e+00   242.182  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 1063/1065
## R-sq.(adj) =  0.992   Deviance explained =   99%
## UBRE = 0.19636  Scale est. = 1         n = 2060

Figura 2.3.2 Predito e observado por sítio

Figura 2.3.3 Diagnostico: avaliação do ajuste pelas preditoras lineares

2.4 Conjunto Completo dos Dados

Dados

The following code was adapted from Gavin Simpson (2016) blog post [acessado em abril de 2020]:
https://fromthebottomoftheheap.net/2016/12/15/simultaneous-interval-revisited/

# new data
df_pred0 <- expand.grid(SiteCode=levels(df_resultados$SiteCode)[1],
                       d.z=seq(min(df_resultados$d.z),max(df_resultados$d.z),length=40),
                       p.z=seq(min(df_resultados$p.z),max(df_resultados$p.z),length=200))
l_mdNref <- list()
l_mdNref[[1]] <- l_nRefEE_mdCheio[["gamm d|Site for each"]]
l_mdNref[[2]] <- l_nRefEI_mdCheio[["gamm d|Site for each"]]
names(l_mdNref) <- c("EE","EI")
# function
f_predict_PostDist <- function(gamm,df=df_pred0,N=1e4){
  #objetos em comum
  beta <- coef(gamm)
  Vb <- vcov(gamm)
  Xp <- predict(gamm,df,type="lpmatrix")
  pred <- predict.gam(gamm,df,se.fit = TRUE)
  se.fit <- pred$se.fit 
  pred.response <- predict.gam(gamm,newdata = df,type = "response",se.fit = TRUE)
  ### estabelece um fator de correção para se.fit
  BUdiff <- MASS::mvrnorm(N,mu=rep(0,nrow(Vb)),Sigma=Vb)
  simDev <- Xp %*% t(BUdiff)
  absDev <- abs(sweep(simDev, 1, se.fit, FUN = "/"))
  masd <- apply(absDev, 2L, max)
  crit <- quantile(masd,prob=0.95,type=8)
  ### simula todos os intervalos de confinça
  ### e então calcula os quantis para cada ponto predito
  sims <- MASS::mvrnorm(N,mu=beta,Sigma=Vb)
  fits <- Xp %*% t(sims)
  q_fits <- t(apply(fits,1,quantile,prob=c(.025,0.5,0.975)))
  ## output
  df_return <- cbind(df,
                     data.frame(fit = pred$fit,
                                se.fit = pred$se.fit,
                                fitResponse = pred.response$fit,
                                se.fitResponse = pred.response$se.fit,
                                avg_fit = apply(fits,1,mean),
                                median_fit = q_fits[,2],
                                q_.025 = q_fits[,1],
                                q_.975 = q_fits[,3])
                     )
  df_return %<>% 
    mutate(upper.fit = fit + (crit*se.fit),
           lower.fit = fit - (crit*se.fit)) %>% 
    select(SiteCode,d.z,p.z,fitResponse,se.fitResponse,fit,se.fit,upper.fit,lower.fit,avg_fit,median_fit,q_.025,q_.975)
  return(df_return)
}
df_plot <- ldply(l_mdNref,f_predict_PostDist,.id="MN")

Figura 2.4.1 Comparação dos métodos de obter os resultados. Na primeira coluna comparação

#
df_plot %<>% mutate(upper.fit_link = fit + (2*se.fit),
                    lower.fit_link = fit + (2*se.fit))
df_plot1 <- as.data.frame(apply(df_plot[,7:16],2,function(X) arm::invlogit(X)*100))
names(df_plot1) <- sapply(names(df_plot1), function(x) paste0(x,"_response"))
df_plot %<>% cbind(.,df_plot1)
#
l_p <- list()
l_p[[1]] <- ggplot(df_plot,aes(x=fitResponse,y=fit_response)) +
  geom_point(alpha=0.3) + labs(title = "mean:invLink(link) ~ response")
l_p[[2]] <- ggplot(df_plot,aes(x=se.fitResponse,y=se.fit_response)) +
  geom_point(alpha=0.3) + labs(title = "se:invLink(link) ~ response")
l_p[[3]] <- ggplot(df_plot,aes(x=fit_response,y=median_fit_response)) +
  geom_point(alpha=0.3) + labs(title = "median posterior ~ link")
l_p[[4]] <- ggplot(df_plot,aes(x=upper.fit_response,y=q_.025_response)) +
  geom_point(alpha=0.3) + labs(title = " 97,5% posterior ~ fit + crit*se.fit")
l_p[[5]] <- ggplot(df_plot,aes(x=lower.fit_response,y=q_.975_response)) +
  geom_point(alpha=0.3) + labs(title = "2,5% posterior ~ fit - crit*se.fit")
# do.call("grid.arrange",c(l_p,ncol=3))
grid.arrange(l_p[[1]],l_p[[2]],l_p[[3]],l_p[[4]],l_p[[5]],
             layout_matrix=rbind(c(1,3,NA),
                                 c(2,4,5)),
             top="response scale")